! ******************************************************
! Eulerian Step and Lagrangian Step
! ******************************************************

module Physics
    use Boundary
    use Communication

    implicit none

contains

    subroutine updateRunTime()
        implicit none
        
        Runtime = Runtime+dt
    end subroutine updateRunTime

    subroutine update_ne_prev()
        implicit none

        ne_prev = ne
    end subroutine update_ne_prev

    subroutine update_ne_temp()
        implicit none

        ne_temp = ( ne_prev + ne )*0.5
    end subroutine update_ne_temp

    subroutine advanceIons()
        implicit none
        type(PseudoParticle), pointer :: p, q
        logical :: IonIsInside, InjectorIsInside

        if(UseInjector)then
            q => InjectorListHead
            p => InjectorListHead%next
            do while(associated(p))
                call pushIons(p)
                call isInInjector(p,InjectorIsInside)
                call deleteIon(InjectorListHead,q)
                p => q%next
            end do
        end if

        q => IonListHead
        p => IonListHead%next
        do while(associated(p))
            call pushIons(p)
            call isInside(p,IonIsInside)
            if(IonIsInside)then
                q => p
                p => p%next
            else
                call deleteIon(IonListHead,q)
                p => q%next
            end if
        end do

        ! Boundaries at Each Direction
        call applyParticleBoundary()  ! Boundary.F90

        ! Exchange Ions on Each Cell
        call exchangeIons()  ! Communication.F90
        if(UseInjector)then
            q => IonListHead
            p => IonListHead%next
            do while(associated(p))
                call isInInjector2(p,InjectorIsInside) ! reload injector
                if(InjectorIsInside) call copyIon(InjectorListHead,p)
                q => p
                p => p%next
            end do
        end if
    end subroutine advanceIons

    subroutine pullbackIons()
        implicit none
        type(PseudoParticle), pointer :: p, q
        logical :: IonIsInside, InjectorIsInside

        if(UseInjector)then
            q => InjectorListHead
            p => InjectorListHead%next
            do while(associated(p))
                call pushIons(p)
                call isInInjector(p,InjectorIsInside)
                call deleteIon(InjectorListHead,q)
                p => q%next
            end do
        end if

        q => IonListHead
        p => IonListHead%next
        do while(associated(p))
            call pushIons(p)
            call isInside(p,IonIsInside)
            if(IonIsInside)then
                q => p
                p => p%next
            else
                call deleteIon(IonListHead,q)
                p => q%next
            end if
        end do

        ! Boundaries at Each Direction
        call applyOpenForParticle()
        call applyPeriodicForParticle()
        call applyReflectForParticle()
        call applyThermalForParticle()  ! Boundary.F90

        ! Exchange Ions on Each Cell
        call exchangeIons()  ! Communication.F90
        if(UseInjector)then
            q => IonListHead
            p => IonListHead%next
            do while(associated(p))
                call isInInjector2(p,InjectorIsInside) ! reload injector
                if(InjectorIsInside) call copyIon(InjectorListHead,p)
                q => p
                p => p%next
            end do
        end if
    end subroutine pullbackIons

    subroutine pushIons(p)
        implicit none
        type(PseudoParticle), pointer :: p
        real(kind=8) :: vn(3), v0(3), vp(3), v(3)
        real(kind=8) :: Efield(3), Bfield(3)
        real(kind=8) :: d(3), s(3)
        real(kind=8) :: eta

        ! Non-relativistic Boris' Method
        v(1) = p%vx
        v(2) = p%vy
        v(3) = p%vz
        call getEBfieldOnParticle(p,Efield,Bfield)
        ! Efield = 0.0d0
        ! Bfield = 0.0d0
        eta = dt*p%Z/(2*p%A*mp/M_unit)

        vn = v+eta*Efield
        d = eta*Bfield
        v0 = vn+cross_product(vn,d)
        s = 2.0d0/(1.0d0+dot_product(d,d))*d
        vp = vn+cross_product(v0,s)
        v = vp+eta*Efield

        p%x = p%x+v(1)*dt
        p%y = p%y+v(2)*dt
        p%vx = v(1) 
        p%vy = v(2) 
        p%vz = v(3)
        
    end subroutine pushIons

    subroutine pullIons(p)
        implicit none
        type(PseudoParticle), pointer :: p
        real(kind=8) :: v(3), M_inv(3,3), M(3,3)
        real(kind=8) :: Efield(3), Bfield(3)
        real(kind=8) :: d(3), s(3)
        real(kind=8) :: eta, ds

        ! Non-relativistic Boris' Method
        v(1) = p%vx
        v(2) = p%vy
        v(3) = p%vz
        p%x = p%x-v(1)*dt
        p%y = p%y-v(2)*dt
        call getEBfieldOnParticle(p,Efield,Bfield)  
        eta = dt*p%Z/(2*p%A*mp/M_unit)
        d = eta*Bfield
        s = 2.0d0/(1.0d0+dot_product(d,d))*d
        ds = dot_product(d,s)
        M = (1-ds)*reshape([1,0,0,0,1,0,0,0,1],[3,3]) &
        + reshape([0.0d0,s(3),-s(2),s(3),0.0d0,-s(1),s(2),-s(1),0.0d0],[3,3]) &
        + reshape([d(1)*s(1),d(1)*s(2),d(1)*s(3),d(2)*s(1),d(2)*s(2),d(2)*s(3),d(3)*s(1),d(3)*s(2),d(3)*s(3)],[3,3])

        M_inv = MAT3_INV(M) ! Setup.F90

        v = v-2*eta*Efield-cross_product(eta*Efield,s)-cross_product(cross_product(eta*Efield,d),s)
        v = matmul(M_inv,v)

        p%vx = v(1) 
        p%vy = v(2) 
        p%vz = v(3)
    end subroutine pullIons

    subroutine getEBfieldOnParticle(p,Efield,Bfield)
        implicit none
        type(PseudoParticle), pointer :: p
        real(kind=8) :: Efield(3), Bfield(3)
        real(kind=8) :: rp(2)
        integer(kind=4) :: Index(2)

        Index(1) = floor((p%x-xMin)/dx)
        Index(2) = floor((p%y-yMin)/dy)
        rp(1) = (xMin-p%x)/dx+Index(1)+1
        rp(2) = (yMin-p%y)/dy+Index(2)+1

        Efield(1) = mapGrid2Particle(rp, &
                    Exg(Index(1):Index(1)+1,Index(2):Index(2)+1))
        Efield(2) = mapGrid2Particle(rp, &
                    Eyg(Index(1):Index(1)+1,Index(2):Index(2)+1))
        Efield(3) = mapGrid2Particle(rp, &
                    Ez(Index(1):Index(1)+1,Index(2):Index(2)+1))
        Bfield(1) = mapGrid2Particle(rp, &
                    (Bxg(Index(1):Index(1)+1,Index(2):Index(2)+1)))
        Bfield(2) = mapGrid2Particle(rp, &
                    (Byg(Index(1):Index(1)+1,Index(2):Index(2)+1)))
        Bfield(3) = mapGrid2Particle(rp, &
                    (Bzg(Index(1):Index(1)+1,Index(2):Index(2)+1)))
    end subroutine getEBfieldOnParticle

    subroutine calculateDistribution()
        implicit none
        type(PseudoParticle), pointer :: p
        real(kind=8) :: rp(2)
        integer(kind=4) :: Index(2),i,j

        ni = 0.0d0
        ne = 0.0d0
        ux = 0.0d0
        uy = 0.0d0
        uz = 0.0d0

        p => IonListHead%next
        do while(associated(p))
            Index(1) = floor((p%x-xMin)/dx)
            Index(2) = floor((p%y-yMin)/dy)
            rp(1) = (xMin-p%x)/dx+Index(1)+1
            rp(2) = (yMin-p%y)/dy+Index(2)+1
            
            ni(Index(1):Index(1)+1,Index(2):Index(2)+1) = &
                ni(Index(1):Index(1)+1,Index(2):Index(2)+1)+ &
                mapParticle2Grid(rp,1.0d0,p%weight,1)
            ne(Index(1):Index(1)+1,Index(2):Index(2)+1) = &
                ne(Index(1):Index(1)+1,Index(2):Index(2)+1)+ &
                mapParticle2Grid(rp,1.0d0,p%weight,p%Z)  ! charge weighted density
            ! Note that here u is ni*u.
            ux(Index(1):Index(1)+1,Index(2):Index(2)+1) = &
                ux(Index(1):Index(1)+1,Index(2):Index(2)+1)+ &
                mapParticle2Grid(rp,p%vx,p%weight,1)
            uy(Index(1):Index(1)+1,Index(2):Index(2)+1) = &
                uy(Index(1):Index(1)+1,Index(2):Index(2)+1)+ &
                mapParticle2Grid(rp,p%vy,p%weight,1)
            uz(Index(1):Index(1)+1,Index(2):Index(2)+1) = &
                uz(Index(1):Index(1)+1,Index(2):Index(2)+1)+ &
                mapParticle2Grid(rp,p%vz,p%weight,1)
            p => p%next
        end do

        ! Exchange ni*u and ni
        call exchangeDistribution()  ! Communication.F90

        ! Boundaries at Each Direction
        call applyGridBoundary()  ! Boundary.F90
        
        ! Calculating Fluid Velocity
        where(ni>0.0d0)
            ux = ux/ni
            uy = uy/ni
            uz = uz/ni
        end where
        
        call smoothGrid(ux)
        call smoothGrid(uy)
        call smoothGrid(uz)
        call smoothGrid(ni)
        call smoothGrid(ne)
    end subroutine calculateDistribution

    subroutine calculateTemperature()
        implicit none

        ! Assuming that the Plasma is an adiabatic ideal gas
        where(ne0>0.0d0) Te = Te0*(ne_prev/ne0)**(Gamma-1)
        call smoothGrid(Te)
    end subroutine calculateTemperature

    function calculatePoissonKernel(ip,jp)
        implicit none
        integer(kind=4) :: i, j, ip, jp, k, Nk
        real(kind=8) :: calculatePoissonKernel(-ncx:ncx,-ncy:ncy)
        real(kind=8) :: LD, KernelCenter, temp, a

        a = sqrt(ne_prev(ip,jp)/Te(ip,jp))
        do j = -ncy,ncy
            do i = -ncx,ncx
                calculatePoissonKernel(i,j) = a**2*besselk0(a*sqrt((i*dx)**2+(j*dy)**2))*dx*dy/2/pi
            end do 
        end do 
        calculatePoissonKernel(0,0) = 0.0d0
        calculatePoissonKernel(0,0) = 1-exp(-a*dx/2.0d0) - sum(calculatePoissonKernel(0,:))
        calculatePoissonKernel = calculatePoissonKernel/sum(calculatePoissonKernel)

    end function calculatePoissonKernel

    function calculateICSKernel(ip,jp,n_ele)
        implicit none
        integer(kind=4) :: i, j, ip, jp, k, Nk
        real(kind=8) :: n_ele(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: calculateICSKernel(-ncx:ncx,-ncy:ncy)
        real(kind=8) :: LD, KernelCenter, temp, a

        a = sqrt(n_ele(ip,jp)/electronTrap)
        
        do j = -ncy,ncy
            do i = -ncx,ncx
                calculateICSKernel(i,j) = a**2*besselk0(a*sqrt((i*dx)**2+(j*dy)**2))*dx*dy/2/pi
            end do 
        end do 
        ! write(*,*) dx,dy
        calculateICSKernel(0,0) = 0.0d0
        calculateICSKernel(0,0) = 1-exp(-a*dx/2.0d0) - sum(calculateICSKernel(0,:))
        calculateICSKernel = calculateICSKernel/sum(calculateICSKernel)

    end function calculateICSKernel

    subroutine calculate_ne_conv()
        implicit none
        integer(kind=4) :: i, j
        real(kind=8) :: lap_Te(-ncx:Nx+ncx,-ncy:Ny+ncy), B_div(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp(-ncx:Nx+ncx,-ncy:Ny+ncy), Bcurl(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: Bcurl1(-ncx:Nx+ncx,-ncy:Ny+ncy), Bcurl2(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: Bcurl3(-ncx:Nx+ncx,-ncy:Ny+ncy), Bcurl4(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: ucrossB(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1,3), B_term(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1,3)
        real(kind=8) :: Temp2(Nx+1+2*ncx,Ny+1+2*ncy), Kernel(-ncx:ncx,-ncy:ncy)

        ! if(UseCAMCL .or. UseField)then
            ! lap_Te = (Te(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy)+Te(1-ncx:Nx+ncx+1,-ncy:Ny+ncy) &
            ! -2*Te(-ncx:Nx+ncx,-ncy:Ny+ncy))/dx/dx &
            ! + (Te(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1)+Te(-ncx:Nx+ncx,1-ncy:Ny+ncy+1) &
            ! -2*Te(-ncx:Nx+ncx,-ncy:Ny+ncy))/dy/dy

        !     ! Temp1 = (Bz(0:Nx,1:Ny+1)-Bz(0:Nx,-1:Ny-1))/2/dy! DyBz-DzBy
        !     ! Temp2 = (Bz(-1:Nx-1,0:Ny)-Bz(1:Nx+1,0:Ny))/2/dx! DzBx-DxBz
        !     ! Temp3 = (By(1:Nx+1,0:Ny)-By(-1:Nx-1,0:Ny))/2/dx-(Bx(0:Nx,1:Ny+1)-Bx(0:Nx,-1:Ny-1))/2/dy! DxBy-DyBx
        !     ! Bcurl = cross_product_arr2(Bx(-ncx:Nx+ncx,-ncy:Ny+ncy) &
        !     !                         ,By(-ncx:Nx+ncx,-ncy:Ny+ncy) &
        !     !                         ,Bz(-ncx:Nx+ncx,-ncy:Ny+ncy) &
        !     !                         ,Temp1,Temp2,Temp3,Nx+1+2*ncx,Ny+1+2*ncy)
        !     Bcurl1 = ((By(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)-By(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy))/2/dx) &
        !             *((By(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)/ne_prev(1-ncx:Nx+ncx+1,-ncy:Ny+ncy) &
        !             -By(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy)/ne_prev(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy))/2/dx) &
        !             +((Bz(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)-Bz(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy))/2/dx) &
        !             *((Bz(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)/ne_prev(1-ncx:Nx+ncx+1,-ncy:Ny+ncy) &
        !             -Bz(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy)/ne_prev(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy))/2/dx) &
        !             +((Bx(-ncx:Nx+ncx,1-ncy:Ny+ncy+1)-Bx(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1))/2/dy) &
        !             *((Bx(-ncx:Nx+ncx,1-ncy:Ny+ncy+1)/ne_prev(-ncx:Nx+ncx,1-ncy:Ny+ncy+1) &
        !             -Bx(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1)/ne_prev(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1))/2/dy) &
        !             +((Bz(-ncx:Nx+ncx,1-ncy:Ny+ncy+1)-Bz(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1))/2/dy) &
        !             *((Bz(-ncx:Nx+ncx,1-ncy:Ny+ncy+1)/ne_prev(-ncx:Nx+ncx,1-ncy:Ny+ncy+1) &
        !             -Bz(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1)/ne_prev(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1))/2/dy)

        !     Temp2 = (Bx*Bx+By*By+Bz*Bz)/ne_prev
        !     Bcurl2 = (Temp2(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy)+Temp2(1-ncx:Nx+ncx+1,-ncy:Ny+ncy) &
        !     -2*Temp2(-ncx:Nx+ncx,-ncy:Ny+ncy))/dx/dx &
        !     + (Temp2(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1)+Temp2(-ncx:Nx+ncx,1-ncy:Ny+ncy+1) &
        !     -2*Temp2(-ncx:Nx+ncx,-ncy:Ny+ncy))/dy/dy

        !     Bcurl3 = ((By(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)/ne_prev(1-ncx:Nx+ncx+1,-ncy:Ny+ncy) &
        !     -By(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy)/ne_prev(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy))/2/dx)  &
        !             *((Bx(-ncx:Nx+ncx,1-ncy:Ny+ncy+1)-Bx(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1))/2/dy) &
        !             +((By(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)-By(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy))/2/dx) &
        !             *((Bx(-ncx:Nx+ncx,1-ncy:Ny+ncy+1)/ne_prev(-ncx:Nx+ncx,1-ncy:Ny+ncy+1) &
        !             -Bx(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1)/ne_prev(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1))/2/dy)

        !     Bcurl4 = (Bx(-ncx:Nx+ncx,-ncy:Ny+ncy)*((Bx(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy) &
        !     +Bx(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)-2*Bx(-ncx:Nx+ncx,-ncy:Ny+ncy))/dx/dx   &
        !     +(By(1-ncx:Nx+ncx+1,1-ncy:Ny+ncy+1)+By(-1-ncx:Nx+ncx-1,-1-ncy:Ny+ncy-1) &
        !     -By(1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy-1)-By(-1-ncx:Nx+ncx-1,1-ncy:Ny+ncy+1))/4/dx/dy) &
        !     +By(-ncx:Nx+ncx,-ncy:Ny+ncy)*((By(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy) &
        !     +By(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)-2*By(-ncx:Nx+ncx,-ncy:Ny+ncy))/dy/dy   &
        !     +(Bx(1-ncx:Nx+ncx+1,1-ncy:Ny+ncy+1)+Bx(-1-ncx:Nx+ncx-1,-1-ncy:Ny+ncy-1) &
        !     -Bx(1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy-1)-Bx(-1-ncx:Nx+ncx-1,1-ncy:Ny+ncy+1))/4/dx/dy) &
        !     )/ne_prev(-ncx:Nx+ncx,-ncy:Ny+ncy)

        !     Bcurl = Bcurl1 + Bcurl2 - Bcurl3 - Bcurl4

        !     ucrossB = cross_product_arr2(ux,uy,uz,Bx,By,Bz,Nx+3+2*ncx,Ny+3+2*ncy)

        !     B_term(:,:,1) = (ne*ucrossB(:,:,1)/ne_prev)
        !     B_term(:,:,2) = (ne*ucrossB(:,:,2)/ne_prev)
        !     B_term(:,:,3) = (ne*ucrossB(:,:,3)/ne_prev)
        !     B_div = ((B_term(1-ncx:Nx+ncx+1,-ncy:Ny+ncy,1)-B_term(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy,1))/2/dx &
        !             +(B_term(-ncx:Nx+ncx,1-ncy:Ny+ncy+1,2)-B_term(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1,2))/2/dy ) &
        !             /ne_prev + Bcurl  !div
        ! else
            lap_Te = 0.0d0
            B_div = 0.0d0
        ! end if
        temp = ne(-ncx:Nx+ncx,-ncy:Ny+ncy) + lap_Te + B_div

        do j=0,Ny 
            do i=0,Nx
                Kernel = calculatePoissonKernel(i,j)
                ne(i,j) = sum(temp(i-ncx:i+ncx,j-ncy:j+ncy)*Kernel)
            end do
        end do
        ! if(myid==1) then
        !     write(*,*) ncx,ncy
        !     do i = -ncx,ncx
        !         write(*,*) Kernel(i,-ncy:ncy)
        !     end do
        ! end if
        ! call smoothGrid(ne)
    end subroutine calculate_ne_conv

    subroutine calculatePressure(n_ele)
        implicit none
        real(kind=8) :: n_ele(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        ! Assuming that the Plasma is an adiabatic ideal gas
        where(ne0>0.0d0) pe = pe0*(n_ele/ne0)**Gamma
        call smoothGrid(pe)
        ! call exchangeGrid(pe)
    end subroutine calculatePressure

    subroutine useElectrostaticOnly()
        implicit none
    ! **** Advancing Ions ******************************
    ! Pushing Every Ions at This Step
        call advanceIons()  ! Physics.F90
    ! **************************************************

    ! **** MPI Communication for Grid ******************
        ! Calculating u and ni at N+1/2 Step
        call calculateDistribution()  ! Physics.F90
    ! **************************************************
        if(UseConv) then
            call calculateTemperature()

            call calculate_ne_conv()
        end if
    ! **** Calculating EOS *****************************
        ! Calculating Pressure by the EOS at N+1 step
        call calculatePressure(ne)  ! Physics.F90
    ! **************************************************
        Ex(1:Nx,0:Ny) = 2*(pe(0:Nx-1,0:Ny)-pe(1:Nx,0:Ny))/(ne(0:Nx-1,0:Ny)+ne(1:Nx,0:Ny))/dx
        Ey(0:Nx,1:Ny) = 2*(pe(0:Nx,0:Ny-1)-pe(0:Nx,1:Ny))/(ne(0:Nx,0:Ny-1)+ne(0:Nx,1:Ny))/dy

        call exchangeGrid(Ex)
        call exchangeGrid(Ey)

        call smoothGrid(Ex)
        call smoothGrid(Ey)
        Exg = mapEfield2Grid(Ex)
        Eyg = mapEfield2Grid(Ey)
    end subroutine useElectrostaticOnly

    subroutine useCAM_CL_Method()
        implicit none
        integer(kind=4) :: i
        real(kind=8) :: utx(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: uty(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: utz(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        call advanceIons()  ! Physics.F90
    ! **************************************************
        do i=1,nt_b
            ! Solving the Faraday Equation
            call solveFaradayEquation(0.5*dt/nt_b,ne_prev)  ! Physics.F90
        end do
    ! **** MPI Communication for Grid ******************
            ! Calculating u and ni at N+1/2 Step
        call calculateDistribution()  ! Physics.F90
    ! **************************************************
        if(UseConv) then
            call calculateTemperature()

            call calculate_ne_conv()
        end if
        call update_ne_temp() ! ne at N+1/2 step
    ! **** Calculating EOS *****************************
        ! Calculating Pressure by the EOS at N+1 step
        call calculatePressure(ne_temp)  ! Physics.F90

    ! **** CAM-CL Method *******************************
        
        ! **** B Field of the Ion-Weibel Instability *******
        ! Solving the ICS Equation
        call calculateIonWeibel(ne_temp)  ! Physics.F90
        ! **************************************************

        ! Calculating the Total Magnetic Field at N+1/2 Step
        call getMagneticField()  ! Physics.F90
        ! Solving the General Ohm Law
        call solveOhmEquation(Bx,By,Bz,ne_temp,ux,uy,uz)  ! Physics.F90
        
        do i=1,nt_b
            ! Solving the Faraday Equation
            call solveFaradayEquation(0.5*dt/nt_b,ne_temp)  ! Physics.F90
        end do
        ! **** B Field of the Ion-Weibel Instability *******
        ! Solving the ICS Equation
        call calculateIonWeibel(ne_temp)  ! Physics.F90
        ! **************************************************

        ! Calculating the Total Magnetic Field at N+1 Step
        call getMagneticField()  ! Physics.F90

        call calculatePressure(ne)  ! Physics.F90
        ! Solving the General Ohm Law
        call solveOhmEquation(Bx,By,Bz,ne,ux,uy,uz)  ! Physics.F90

        ! Solving the Lorentz Equation
        ! Solving the General Ohm Law
        utx = ux
        uty = uy
        utz = uz
        call solveLorentzEquation() 

        call solveOhmEquation(Bx,By,Bz,ne,ux,uy,uz)
        ux = utx
        uy = uty
        uz = utz
    ! **************************************************
        call update_ne_prev()

        call exchangeGrid(Bx)
        call exchangeGrid(By)
        call exchangeGrid(Bz)

        call smoothGrid(Bx)
        call smoothGrid(By)
        call smoothGrid(Bz)

        call exchangeGrid(Ex)
        call exchangeGrid(Ey)
        call exchangeGrid(Ez)

        call smoothGrid(Ex)
        call smoothGrid(Ey)
        call smoothGrid(Ez)

    end subroutine useCAM_CL_Method

    subroutine calculateIonWeibel(n_ele)
        implicit none
        real(kind=8) :: n_ele(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        if(UseICS)then
            call solveICSEquation_conv(n_ele)
            ! call solveICSEquation_SI(n_ele)
        else
            Bwx = 0.0d0
            Bwy = 0.0d0
            Bwz = 0.0d0
        end if
    end subroutine calculateIonWeibel

    subroutine solveICSEquation_conv(n_ele)
        implicit none
        integer(kind=4) :: i,j
        real(kind=8) :: n_ele(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Temp1(-ncx:Nx+ncx,-ncy:Ny+ncy), Temp2(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: Temp3(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: Tempx(0:Nx,0:Ny), Tempy(0:Nx,0:Ny), Tempz(0:Nx,0:Ny)
        real(kind=8) :: temp_bx(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_by(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_bz(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: Kernel(-ncx:ncx,-ncy:ncy)

        temp_bx = mapBfield2Grid(Bfx)
        temp_by = mapBfield2Grid(Bfy)
        temp_bz = mapBfield2Grid(Bfz)

        Temp1 = electronTrap*(uz(-ncx:Nx+ncx,1-ncy:Ny+ncy+1)-uz(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1))/2/dy + temp_bx! Dyuz-Dzuy
        Temp2 = electronTrap*(uz(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy)-uz(1-ncx:Nx+ncx+1,-ncy:Ny+ncy))/2/dx + temp_by! Dzux-Dxuz
        Temp3 = electronTrap*((uy(1-ncx:Nx+ncx+1,-ncy:Ny+ncy)-uy(-1-ncx:Nx+ncx-1,-ncy:Ny+ncy))/2/dx &
        - (ux(-ncx:Nx+ncx,1-ncy:Ny+ncy+1)-ux(-ncx:Nx+ncx,-1-ncy:Ny+ncy-1))/2/dy) + temp_bz! Dxuy-Dyux

        do j=0,Ny 
            do i=0,Nx
                Kernel = calculateICSKernel(i,j,n_ele)
                Tempx(i,j) = sum(Temp1(i-ncx:i+ncx,j-ncy:j+ncy)*Kernel)
                Tempy(i,j) = sum(Temp2(i-ncx:i+ncx,j-ncy:j+ncy)*Kernel)
                Tempz(i,j) = sum(Temp3(i-ncx:i+ncx,j-ncy:j+ncy)*Kernel)
            end do
        end do
        ! if(myid==1) then
        !     ! write(*,*) ncx,ncy,Kernel(0,0)
        !     do i = -ncx,ncx
        !         write(*,*) Kernel(i,-ncy:ncy)
        !     end do
        ! end if
        Bwx(0:Nx,1:Ny) = (Tempx(0:Nx,1:Ny)+Tempx(0:Nx,0:Ny-1))/2
        Bwy(1:Nx,0:Ny) = (Tempy(1:Nx,0:Ny)+Tempy(0:Nx-1,0:Ny))/2
        Bwz(1:Nx,1:Ny) = (Tempz(1:Nx,1:Ny)+Tempz(0:Nx-1,1:Ny)+Tempz(1:Nx,0:Ny-1)+Tempz(0:Nx-1,0:Ny-1))/4
        
        ! Bwx = 0.0d0
        ! Bwy = 0.0d0
        ! Bwz = 0.0d0
        call exchangeGrid(Bwx)
        call exchangeGrid(Bwy)
        call exchangeGrid(Bwz)
    end subroutine solveICSEquation_conv

    subroutine solveICSEquation_SI(n_ele)
        implicit none
        integer(kind=4) :: i
        real(kind=8) :: n_ele(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: invex(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), invey(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: invez(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: temp_nex(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), temp_ney(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: temp_nez(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: temp_bx(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), temp_by(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: temp_bz(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: Temp1(0:Nx,0:Ny), Temp2(0:Nx,0:Ny), Temp3(0:Nx,0:Ny)

        invex = 2.0/dx/dx+2.0/dy/dy+mapGrid2BxGrid(n_ele)/electronTrap ! 2D case
        invey = 2.0/dx/dx+2.0/dy/dy+mapGrid2ByGrid(n_ele)/electronTrap
        invez = 2.0/dx/dx+2.0/dy/dy+mapGrid2BzGrid(n_ele)/electronTrap
        !in first calculation need ne pe u at N+1/2 step

        Temp1 = (uz(0:Nx,1:Ny+1)-uz(0:Nx,-1:Ny-1))/2/dy! Dyuz-Dzuy
        Temp2 = (uz(-1:Nx-1,0:Ny)-uz(1:Nx+1,0:Ny))/2/dx! Dzux-Dxuz
        Temp3 = (uy(1:Nx+1,0:Ny)-uy(-1:Nx-1,0:Ny))/2/dx - (ux(0:Nx,1:Ny+1)-ux(0:Nx,-1:Ny-1))/2/dy! Dxuy-Dyux
        ! Simplified Version
        ! Bwx(0:Nx,1:Ny) = (Temp1(0:Nx,1:Ny)+Temp1(0:Nx,0:Ny-1))/2
        ! Bwy(1:Nx,0:Ny) = (Temp2(1:Nx,0:Ny)+Temp2(0:Nx-1,0:Ny))/2
        ! Bwz(1:Nx,1:Ny) = (Temp3(1:Nx,1:Ny)+Temp3(0:Nx-1,1:Ny)+Temp3(1:Nx,0:Ny-1)+Temp3(0:Nx-1,0:Ny-1))/4
        ! Simple iteration
        do i=1,8
            temp_bx = getBxGrid(Bwx)
            temp_by = getByGrid(Bwy)
            temp_bz = getBzGrid(Bwz)
            
            Bwx(0:Nx,1:Ny) = ((temp_bx(1:Nx+1,1:Ny)+temp_bx(-1:Nx-1,1:Ny))/dx/dx &
                            +(temp_bx(0:Nx,2:Ny+1)+temp_bx(0:Nx,0:Ny-1))/dy/dy &
                            +n_ele(0:Nx,1:Ny)*Temp1(0:Nx,1:Ny)/2  &
                            +n_ele(0:Nx,0:Ny-1)*Temp1(0:Nx,0:Ny-1)/2)/invex(0:Nx,1:Ny)
            Bwy(1:Nx,0:Ny) = ((temp_by(2:Nx+1,0:Ny)+temp_by(0:Nx-1,0:Ny))/dx/dx &
                            +(temp_by(1:Nx,1:Ny+1)+temp_by(1:Nx,-1:Ny-1))/dy/dy &
                            +n_ele(1:Nx,0:Ny)*Temp2(1:Nx,0:Ny)/2  &
                            +n_ele(1:Nx,0:Ny)*Temp2(0:Nx-1,0:Ny)/2)/invey(1:Nx,0:Ny)
            Bwz(1:Nx,1:Ny) = ((temp_bz(2:Nx+1,1:Ny)+temp_bz(0:Nx-1,1:Ny))/dx/dx &
                            +(temp_bz(1:Nx,2:Ny+1)+temp_bz(1:Nx,0:Ny-1))/dy/dy &
                            +n_ele(1:Nx,1:Ny)*Temp3(1:Nx,1:Ny)/4  &
                            +n_ele(0:Nx-1,1:Ny)*Temp3(0:Nx-1,1:Ny)/4  &
                            +n_ele(1:Nx,0:Ny-1)*Temp3(1:Nx,0:Ny-1)/4  &
                            +n_ele(0:Nx-1,0:Ny-1)*Temp3(0:Nx-1,0:Ny-1)/4)/invez(1:Nx,1:Ny)
            call exchangeGrid(Bwx)
            call exchangeGrid(Bwy)
            call exchangeGrid(Bwz)
        end do

        ! Bwx = 0.0d0
        ! Bwy = 0.0d0
        ! Bwz = 0.0d0
        call exchangeGrid(Bwx)
        call exchangeGrid(Bwy)
        call exchangeGrid(Bwz)
    end subroutine solveICSEquation_SI

    subroutine solveFaradayEquation(Dt,n_ele)
        implicit none
        real(kind=8) :: Hall1(Nx+1,Ny), Hall2(Nx,Ny+1), Hall3(Nx,Ny)
        real(kind=8) :: Temp1(Nx+1,Ny+1), Temp2(Nx+1,Ny+1), Temp3(Nx+1,Ny+1)
        real(kind=8) :: part1(0:Nx,0:Ny,3), vez(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: vex(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), vey(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1,3)
        real(kind=8) :: Tenser(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1,3)
        real(kind=8) :: Dt, n_ele(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: pe_b(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: NonAb1(Nx+1,Ny), NonAb2(Nx,Ny+1), NonAb3(Nx,Ny), ne_b(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: temp_bx(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_by(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_bz(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: Tenser1(0:Nx,0:Ny)
        real(kind=8) :: Tenser2(0:Nx,0:Ny)
        real(kind=8) :: Tenser3(0:Nx,0:Ny)
        real(kind=8) :: Dinvne1(0:Nx,0:Ny)
        real(kind=8) :: Dinvne2(0:Nx,0:Ny)
        real(kind=8) :: Dinvne3(0:Nx,0:Ny)
        real(kind=8) :: Ex_b(1:Nx,0:Ny), Ey_b(0:Nx,1:Ny), Ez_b(0:Nx,0:Ny)
        real(kind=8) :: bgx(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), bgy(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: bgz(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), V1(0:Nx,0:Ny), V2(0:Nx,0:Ny)
        real(kind=8) :: udotDB1(Nx+1,Ny), udotDB2(Nx,Ny+1), udotDB3(Nx,Ny)
        real(kind=8) :: curl(Nx+1,Ny+1,3)

        temp_bx = mapBfield2Grid(Bfx)
        temp_by = mapBfield2Grid(Bfy)
        temp_bz = mapBfield2Grid(Bfz)
        bgx = getBxGrid(Bx)
        bgy = getByGrid(By)
        bgz = getBzGrid(Bz)

        vex(0:Nx,0:Ny) = Z_avg*ni(0:Nx,0:Ny)/n_ele(0:Nx,0:Ny)*ux(0:Nx,0:Ny) &
        -((bgz(0:Nx,1:Ny+1)-bgz(0:Nx,0:Ny)+bgz(1:Nx+1,1:Ny+1)-bgz(1:Nx+1,0:Ny))/2/dy)/n_ele(0:Nx,0:Ny)
        vey(0:Nx,0:Ny) = Z_avg*ni(0:Nx,0:Ny)/n_ele(0:Nx,0:Ny)*uy(0:Nx,0:Ny) &
        -((bgz(0:Nx,0:Ny)-bgz(1:Nx+1,0:Ny)+bgz(0:Nx,1:Ny+1)-bgz(1:Nx+1,1:Ny+1))/2/dx)/n_ele(0:Nx,0:Ny)
        vez(0:Nx,0:Ny) = Z_avg*ni(0:Nx,0:Ny)/n_ele(0:Nx,0:Ny)*uz(0:Nx,0:Ny) &
        -((bgy(1:Nx+1,0:Ny)-bgy(0:Nx,0:Ny))/dx - (bgx(0:Nx,1:Ny+1)-bgx(0:Nx,0:Ny))/dy)/n_ele(0:Nx,0:Ny)
        ! vex(0:Nx,0:Ny) = ux(0:Nx,0:Ny) &
        ! -((bgz(0:Nx,1:Ny+1)-bgz(0:Nx,0:Ny)+bgz(1:Nx+1,1:Ny+1)-bgz(1:Nx+1,0:Ny))/2/dy)/n_ele(0:Nx,0:Ny)
        ! vey(0:Nx,0:Ny) = uy(0:Nx,0:Ny) &
        ! -((bgz(0:Nx,0:Ny)-bgz(1:Nx+1,0:Ny)+bgz(0:Nx,1:Ny+1)-bgz(1:Nx+1,1:Ny+1))/2/dx)/n_ele(0:Nx,0:Ny)
        ! vez(0:Nx,0:Ny) = uz(0:Nx,0:Ny) &
        ! -((bgy(1:Nx+1,0:Ny)-bgy(0:Nx,0:Ny))/dx - (bgx(0:Nx,1:Ny+1)-bgx(0:Nx,0:Ny))/dy)/n_ele(0:Nx,0:Ny)
        
        Temp(0:Nx,0:Ny,:) = cross_product_arr2(vex(0:Nx,0:Ny),vey(0:Nx,0:Ny), &
        vez(0:Nx,0:Ny),temp_bx(0:Nx,0:Ny),temp_by(0:Nx,0:Ny),temp_bz(0:Nx,0:Ny),Nx+1,Ny+1)

        call exchangeGrid(Temp(:,:,1))
        call exchangeGrid(Temp(:,:,2))
        call exchangeGrid(Temp(:,:,3))
        call smoothGridOnly(Temp(:,:,1))
        call smoothGridOnly(Temp(:,:,2))
        call smoothGridOnly(Temp(:,:,3))
        call exchangeGrid(Temp(:,:,1))
        call exchangeGrid(Temp(:,:,2))
        call exchangeGrid(Temp(:,:,3))

        ! Ex_b(1:Nx,0:Ny) = - (Temp(0:Nx-1,0:Ny,1)+Temp(1:Nx,0:Ny,1))/2 &
        !                   + 2*(B2(1:Nx,0:Ny)-B2(0:Nx-1,0:Ny))/(n_ele(1:Nx,0:Ny)+n_ele(0:Nx-1,0:Ny))/2/dx
        ! Ey_b(0:Nx,1:Ny) = - (Temp(0:Nx,0:Ny-1,2)+Temp(0:Nx,1:Ny,2))/2 &
        !                   + 2*(B2(0:Nx,1:Ny)-B2(0:Nx,0:Ny-1))/(n_ele(0:Nx,1:Ny)+n_ele(0:Nx,0:Ny-1))/2/dy
        ! Ez_b(0:Nx,0:Ny) = - Temp(0:Nx,0:Ny,3)

        Ex_b(1:Nx,0:Ny) = - (Temp(0:Nx-1,0:Ny,1)+Temp(1:Nx,0:Ny,1))/2 
        Ey_b(0:Nx,1:Ny) = - (Temp(0:Nx,0:Ny-1,2)+Temp(0:Nx,1:Ny,2))/2 
        Ez_b(0:Nx,0:Ny) = - Temp(0:Nx,0:Ny,3)

        NonAb1 = 0.0d0  ! DyneDzpe-DzneDype
        NonAb2 = 0.0d0  ! DzneDxpe-DxneDype
        NonAb3 = ((ne_b(2:Nx+1,1:Ny)-ne_b(0:Nx-1,1:Ny))*(pe_b(1:Nx,2:Ny+1)-pe_b(1:Nx,0:Ny-1)) &
        -(pe_b(2:Nx+1,1:Ny)-pe_b(0:Nx-1,1:Ny))*(ne_b(1:Nx,2:Ny+1)-ne_b(1:Nx,0:Ny-1)))/4/dx/dy/ne_b(1:Nx,1:Ny)**2
        NonAb3 = 0.0d0 ! FD method cannot automaticly be 0 in adiabatic condition, so we have to restrict it to 0

        Temp1 = (Ez_b(0:Nx,1:Ny)-Ez_b(0:Nx,0:Ny-1))/dy! DyEz-DzEy
        Temp2 = (Ez_b(0:Nx-1,0:Ny)-Ez_b(1:Nx,0:Ny))/dx! DzEx-DxEz
        Temp3 = (Ey_b(1:Nx,1:Ny)-Ey_b(0:Nx-1,1:Ny))/dx &
                - (Ex_b(1:Nx,1:Ny)-Ex_b(1:Nx,0:Ny-1))/dy   ! DxEy-DyEx

        Bfx(0:Nx,1:Ny) = Bfx(0:Nx,1:Ny) - Dt*(Temp1+NonAb1)
        Bfy(1:Nx,0:Ny) = Bfy(1:Nx,0:Ny) - Dt*(Temp2+NonAb2)
        Bfz(1:Nx,1:Ny) = Bfz(1:Nx,1:Ny) - Dt*(Temp3+NonAb3)

        call exchangeGrid(Bfx)
        call exchangeGrid(Bfy)
        call exchangeGrid(Bfz)
    end subroutine solveFaradayEquation

    subroutine getMagneticField()
        implicit none
        if(UseICS)then
            Bx = Bwx
            By = Bwy
            Bz = Bwz
        else
            Bx = Bfx + Bwx
            By = Bfy + Bwy
            Bz = Bfz + Bwz
        end if

        call exchangeGrid(Bx)
        call exchangeGrid(By)
        call exchangeGrid(Bz)

        Bxg = mapBfield2Grid(Bx)
        Byg = mapBfield2Grid(By)
        Bzg = mapBfield2Grid(Bz)
    end subroutine getMagneticField

    subroutine solveOhmEquation(Bfield1,Bfield2,Bfield3,n_ele,velx,vely,velz)
        implicit none
        real(kind=8) :: Bfield1(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bfield2(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bfield3(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: temp_bx(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_by(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_bz(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: velx(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: vely(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: velz(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: n_ele(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1,3), vez(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: vex(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1), vey(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: bgx(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), bgy(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: bgz(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), gradne(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1,3)
        real(kind=8) :: Temp1(Nx+1,Ny+1), Temp2(Nx+1,Ny+1), Temp3(Nx+1,Ny+1)
        real(kind=8) :: curl(Nx+1,Ny+1,3), J(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1,3)
    
        temp_bx = mapBfield2Grid(Bfield1)
        temp_by = mapBfield2Grid(Bfield2)
        temp_bz = mapBfield2Grid(Bfield3)
        bgx = getBxGrid(Bfield1)
        bgy = getByGrid(Bfield2)
        bgz = getBzGrid(Bfield3)

        vex(0:Nx,0:Ny) = Z_avg*ni(0:Nx,0:Ny)/n_ele(0:Nx,0:Ny)*ux(0:Nx,0:Ny) &
        -((bgz(0:Nx,1:Ny+1)-bgz(0:Nx,0:Ny)+bgz(1:Nx+1,1:Ny+1)-bgz(1:Nx+1,0:Ny))/2/dy)/n_ele(0:Nx,0:Ny)
        vey(0:Nx,0:Ny) = Z_avg*ni(0:Nx,0:Ny)/n_ele(0:Nx,0:Ny)*uy(0:Nx,0:Ny) &
        -((bgz(0:Nx,0:Ny)-bgz(1:Nx+1,0:Ny)+bgz(0:Nx,1:Ny+1)-bgz(1:Nx+1,1:Ny+1))/2/dx)/n_ele(0:Nx,0:Ny)
        vez(0:Nx,0:Ny) = Z_avg*ni(0:Nx,0:Ny)/n_ele(0:Nx,0:Ny)*uz(0:Nx,0:Ny) &
        -((bgy(1:Nx+1,0:Ny)-bgy(0:Nx,0:Ny))/dx - (bgx(0:Nx,1:Ny+1)-bgx(0:Nx,0:Ny))/dy)/n_ele(0:Nx,0:Ny)
        ! vex(0:Nx,0:Ny) = ux(0:Nx,0:Ny) &
        ! -((bgz(0:Nx,1:Ny+1)-bgz(0:Nx,0:Ny)+bgz(1:Nx+1,1:Ny+1)-bgz(1:Nx+1,0:Ny))/2/dy)/n_ele(0:Nx,0:Ny)
        ! vey(0:Nx,0:Ny) = uy(0:Nx,0:Ny) &
        ! -((bgz(0:Nx,0:Ny)-bgz(1:Nx+1,0:Ny)+bgz(0:Nx,1:Ny+1)-bgz(1:Nx+1,1:Ny+1))/2/dx)/n_ele(0:Nx,0:Ny)
        ! vez(0:Nx,0:Ny) = uz(0:Nx,0:Ny) &
        ! -((bgy(1:Nx+1,0:Ny)-bgy(0:Nx,0:Ny))/dx - (bgx(0:Nx,1:Ny+1)-bgx(0:Nx,0:Ny))/dy)/n_ele(0:Nx,0:Ny)

        ! call exchangeGrid(vex)
        ! call exchangeGrid(vey)
        ! call exchangeGrid(vez)
        ! call smoothGridOnly(vex)
        ! call smoothGridOnly(vey)
        ! call smoothGridOnly(vez)
        ! call exchangeGrid(vex)
        ! call exchangeGrid(vey)
        ! call exchangeGrid(vez)
        
        Temp(0:Nx,0:Ny,:) = cross_product_arr2(vex(0:Nx,0:Ny),vey(0:Nx,0:Ny), &
        vez(0:Nx,0:Ny),temp_bx(0:Nx,0:Ny),temp_by(0:Nx,0:Ny),temp_bz(0:Nx,0:Ny),Nx+1,Ny+1)
        
        ! call exchangeGrid(Temp(:,:,1))
        ! call exchangeGrid(Temp(:,:,2))
        ! call exchangeGrid(Temp(:,:,3))
        ! call smoothGridOnly(Temp(:,:,1))
        ! call smoothGridOnly(Temp(:,:,2))
        ! call smoothGridOnly(Temp(:,:,3))
        ! call exchangeGrid(Temp(:,:,1))
        ! call exchangeGrid(Temp(:,:,2))
        ! call exchangeGrid(Temp(:,:,3))
        
        gradne(0:Nx,0:Ny,1) = (pe(-1:Nx-1,0:Ny)-pe(1:Nx+1,0:Ny))/(n_ele(0:Nx,0:Ny))/2/dx
        gradne(0:Nx,0:Ny,2) = (pe(0:Nx,-1:Ny-1)-pe(0:Nx,1:Ny+1))/(n_ele(0:Nx,0:Ny))/2/dy
        gradne(0:Nx,0:Ny,3) = 0.0d0
        
        Ex(1:Nx,0:Ny) = (gradne(0:Nx-1,0:Ny,1)+gradne(1:Nx,0:Ny,1))/2 - (Temp(0:Nx-1,0:Ny,1)+Temp(1:Nx,0:Ny,1))/2 
        Ey(0:Nx,1:Ny) = (gradne(0:Nx,0:Ny-1,2)+gradne(0:Nx,1:Ny,2))/2 - (Temp(0:Nx,0:Ny-1,2)+Temp(0:Nx,1:Ny,2))/2 
        Ez(0:Nx,0:Ny) = gradne(0:Nx,0:Ny,3) - Temp(0:Nx,0:Ny,3)
    
        call exchangeGrid(Ex)
        call exchangeGrid(Ey)
        call exchangeGrid(Ez)
        Exg = mapEfield2Grid(Ex)
        Eyg = mapEfield2Grid(Ey)
    end subroutine solveOhmEquation

    subroutine solveLorentzEquation()   ! uncompleted
        implicit none
        real(kind=8) :: Temp0(-ncx:Nx+ncx,-ncy:Ny+ncy,3)
        real(kind=8) :: eta
        real(kind=8) :: Temp1(Nx+1,Ny+1), Temp2(Nx+1,Ny+1), Temp3(Nx+1,Ny+1)
        Temp0 = cross_product_arr2(ux(-ncx:Nx+ncx,-ncy:Ny+ncy),uy(-ncx:Nx+ncx,-ncy:Ny+ncy), &
        uz(-ncx:Nx+ncx,-ncy:Ny+ncy),Bxg,Byg,Bzg,Nx+2*ncx+1,Ny+2*ncy+1)

        eta = dt*Z_avg/(2*A_avg*mp/M_unit)  ! step half
        ux(0:Nx,0:Ny) = ux(0:Nx,0:Ny) + eta*(Exg(0:Nx,0:Ny) + Temp0(0:Nx,0:Ny,1))
        uy(0:Nx,0:Ny) = uy(0:Nx,0:Ny) + eta*(Eyg(0:Nx,0:Ny) + Temp0(0:Nx,0:Ny,2))
        uz(0:Nx,0:Ny) = uz(0:Nx,0:Ny) + eta*(Ez(0:Nx,0:Ny) + Temp0(0:Nx,0:Ny,3))

        call exchangeGrid(ux)
        call exchangeGrid(uy)
        call exchangeGrid(uz)
    end subroutine solveLorentzEquation

    subroutine sovleLastStepofCAMCL()
        implicit none
        real(kind=8) :: Temp0(-ncx:Nx+ncx,-ncy:Ny+ncy,3)
        real(kind=8) :: utx(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: uty(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: utz(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: temp_bx(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_by(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_bz(-ncx:Nx+ncx,-ncy:Ny+ncy)
        real(kind=8) :: eta
        real(kind=8) :: Temp1(Nx+1,Ny+1), Temp2(Nx+1,Ny+1), Temp3(Nx+1,Ny+1)
        
        Temp0 = cross_product_arr2(ux(-ncx:Nx+ncx,-ncy:Ny+ncy),uy(-ncx:Nx+ncx,-ncy:Ny+ncy), &
        uz(-ncx:Nx+ncx,-ncy:Ny+ncy),Bxg,Byg,Bzg,Nx+2*ncx+1,Ny+2*ncy+1)

        utx = ux
        uty = uy
        utz = uz
        eta = dt*Z_avg/(2*A_avg*mp/M_unit)  ! step half
        ux(0:Nx,0:Ny) = ux(0:Nx,0:Ny) + eta*(Exg(0:Nx,0:Ny) + Temp0(0:Nx,0:Ny,1))
        uy(0:Nx,0:Ny) = uy(0:Nx,0:Ny) + eta*(Eyg(0:Nx,0:Ny) + Temp0(0:Nx,0:Ny,2))
        uz(0:Nx,0:Ny) = uz(0:Nx,0:Ny) + eta*(Ez(0:Nx,0:Ny) + Temp0(0:Nx,0:Ny,3))
        call exchangeGrid(ux)
        call exchangeGrid(uy)
        call exchangeGrid(uz)
    
        call solveOhmEquation(Bx,By,Bz,ne,ux,uy,uz)
        ux = utx
        uy = uty
        uz = utz
    end subroutine sovleLastStepofCAMCL

    subroutine useCN_PPC_Method()
        implicit none
        real(kind=8) :: Ex_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Ey_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Ez_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bx_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: By_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bz_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        Ex_temp = 0.0d0
        Ey_temp = 0.0d0
        Ez_temp = 0.0d0
        Bx_temp = 0.0d0
        By_temp = 0.0d0
        Bz_temp = 0.0d0
        Ex_temp = Ex    ! store N step value
        Ey_temp = Ey
        Ez_temp = Ez
        Bx_temp = Bx
        By_temp = By
        Bz_temp = Bz

        call calculateEBfield_in_CN()
        call getEBfield_at_Half_step(Ex_temp,Ey_temp,Ez_temp,Bx_temp,By_temp,Bz_temp)

        call advanceIons()  

        call calculateDistribution()
        if(UseConv) then
            call calculateTemperature()

            call calculate_ne_conv()
        end if
        call calculatePressure(ne)

        call pullbackIons()

        call calculateEBfield_in_CN()
        call getEBfield_at_Half_step(Ex_temp,Ey_temp,Ez_temp,Bx_temp,By_temp,Bz_temp)

        call advanceIons()  

        call calculateDistribution()
        if(UseConv) then
            call calculateTemperature()

            call calculate_ne_conv()
        end if

        call calculatePressure(ne)
        call solveLastStepofCNPPC()
        if(UseConv) then
            call update_ne_prev()
        end if

        call exchangeGrid(Bx)
        call exchangeGrid(By)
        call exchangeGrid(Bz)

        call smoothGrid(Bx)
        call smoothGrid(By)
        call smoothGrid(Bz)

        call exchangeGrid(Ex)
        call exchangeGrid(Ey)
        call exchangeGrid(Ez)

        call smoothGrid(Ex)
        call smoothGrid(Ey)
        call smoothGrid(Ez)
    end subroutine

    subroutine calculateEBfield_in_CN()
        implicit none
        real(kind=8) :: Bfx_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bfy_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bfz_temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
    
        call update_ne_temp()
        call calculateIonWeibel(ne_temp)  ! Physics.F90
        ! N+1 step
        ! Bfx_temp = Bfx
        ! Bfy_temp = Bfy
        ! Bfz_temp = Bfz
        call solveFaradayEquation(dt,ne_temp)  ! Physics.F90
        
        call getMagneticField()  ! Physics.F90
        ! Bfx = Bfx_temp
        ! Bfy = Bfy_temp
        ! Bfz = Bfz_temp

        call solveOhmEquation(Bx,By,Bz,ne,ux,uy,uz)  ! Physics.F90
        
    end subroutine calculateEBfield_in_CN

    subroutine solveLastStepofCNPPC()
        implicit none
    
        call update_ne_temp()
        call calculateIonWeibel(ne_temp)  ! Physics.F90
        ! N+1 step
        call solveFaradayEquation(dt,ne_temp)  ! Physics.F90
        
        call getMagneticField()  ! Physics.F90

        call solveOhmEquation(Bx,By,Bz,ne,ux,uy,uz)  ! Physics.F90

        ! call solveGaussEquation()

        ! if(UseConv) then
        !     call calculateTemperature()

        !     call calculate_ne_conv()
        ! end if
        ! call calculatePressure(ne)

        ! call solveOhmEquation(Bx,By,Bz,ne,ux,uy,uz)  ! Physics.F90
        
    end subroutine solveLastStepofCNPPC

    subroutine getEBfield_at_Half_step(Eprev_x,Eprev_y,Eprev_z,Bprev_x,Bprev_y,Bprev_z)
        implicit none
        real(kind=8) :: Eprev_x(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Eprev_y(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Eprev_z(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bprev_x(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bprev_y(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bprev_z(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
    
        Ex = 0.5*(Eprev_x + Ex)
        Ey = 0.5*(Eprev_y + Ey)
        Ez = 0.5*(Eprev_z + Ez)
        Bx = 0.5*(Bprev_x + Bx)
        By = 0.5*(Bprev_y + By)
        Bz = 0.5*(Bprev_z + Bz)

        Exg = mapEfield2Grid(Ex)
        Eyg = mapEfield2Grid(Ey)

        Bxg = mapBfield2Grid(Bx)
        Byg = mapBfield2Grid(By)
        Bzg = mapBfield2Grid(Bz)
    end subroutine getEBfield_at_Half_step

    subroutine solveGaussEquation()
        implicit none

        ne_prev(0:Nx,0:Ny) = Z_avg*ni(0:Nx,0:Ny) - &
                             (Ex(1:Nx+1,0:Ny)-Ex(0:Nx,0:Ny))/dx + &
                             (Ey(0:Nx,1:Ny+1)-Ey(0:Nx,0:Ny))/dy
    end subroutine solveGaussEquation
end module Physics

